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ABSTRACT 

The gravitational clustering of collisionless particles in an expanding universe is mod- 
elled using some simple physical ideas. I show that it is indeed possible to understand 
the nonlinear clustering in terms of three well defined regimes: (1) linear regime (2) 
quasilinear regime which is dominated by scale-invariant radial infall and (3) nonlinear 
regime dominated by nonradial motions and mergers. Modelling each of these regimes 
separately I show how the nonlinear two point correlation function can be related to the 
linear correlation function in heirarchical models. This analysis leads to results which 
are in good agreement with numerical simulations thereby providing an explanation for 
numerical results. The ideas presented here will also serve as a powerful anlytical tool 
to investigate nonlinear clustering in different models. Several implications of the result 
are discussed. 
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The driving force behind the formation of large scale 
structures in the universe is the gravitational field produced 
by density fluctuations. Overdense regions accrete matter at 
the expense of underdense regions allowing inhomogeneities 
in the universe to grow. Observations suggest that the ma- 
terial content of the universe is dominated by dark matter, 
likely to be made of collisionless elementary particles. In that 
case, the gravitational force is mainly due to these particles 
and, to first approximation, we can ignore the complications 
arising from baryonic physics. The evolution of density per- 
turbati ons is then governed purely by the gravitational force. 



on the well known behaviour of a spherically symmetric over- 
dense region in the universe. In the behaviour of such a re- 
gion, one can identify three different regimes of interest: (1) 
In the early stages of the evolution, when the density con- 
trast is small, the evolution is described by linear theory. (2) 
Each of the spherical shells with an initial radius Xi can be 
parametersed by a mass contained inside the shell, M{xi), 
and the energy, E(xi) for the particular shell. Each shell will 
expand to a maximum radius x max oc M/\E\ and then turn 
around and collapse. Such a spherical collaps e and result 



ing evolution al lows a self similar description ( Filmore and 
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ble to study their evolution using linear theory. But once the 
density contrast becomes comparable to unity, linear pertur- 
bation theory breaks down and one must use N-body sim- 
ulations to study the growth of perturbations. While these 
simulations are of some value in making concrete predic- 
tions for specific models, they do not provide clear physical 
insight into the process of non-linear gravitational dynamics. 
To obtain such an insight into this complex problem, it is 
necessary to model the gravitational clustering of collision- 
less particles using simple physical concepts. I shall develop 
one such model in this paper, which - in spite of extreme 
simplicity - reproduces the simulation results for hierarchical 
models fairly accurately. Further, this model also provides 
insight into the clustering process and can be modified to 
take into account more complicated situations. 

The paradigm for understanding the clustering is based 



Goldreich 1984 ) in which each shell act s as though it has a n 
effective radius proportional to x ma x ( Bertschinger 1985). 
This will be the quasilinear phase. (3) The spherical evolu- 
tion will break down during the later stages due to several 
reasons. First of all, non radial motions will arise due to am- 
plification of deviations from spherical symmetry. Secondly, 
the existence of substructure will influence the evolution in 
a non-spherically symmetric way. Finally, in the real uni- 
verse, there will be merging of such clusters [each of which 
could have been centres of spherical overdense regions in the 
begining] which will again destroy the spherical symmetry. 
This will be the nonlinear phase. 

The description given above is sufficiently vague and 
sufficiently well known that one may suspect it can not lead 
to any insight into the problem. In particular, real universe 
is hardly spherical. I will show that it is, however, possible 
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to model the above process in a manner which allows direct 
generalisation to the real universe. 

To do this we will begin by studying the evolution of 
system starting from a gaussian initial fluctuations with 
an initial power spectrum, Pi„(k). The fourier transform of 
the power spectrum defines the correlation function £ (a, x) 
where a oc t 2 ^ 3 is the expansion factor in a universe with 
O = 1. It is more convenient for our purpose to work with 
the average correlation function inside a sphere of radius x, 
defined by 



€(a,x) = Z(a,y)y dy 



(1) 



In the linear regime we have £l(o,,x) oc a £i n (ai,x). In the 
quasilinear and nonlinear regimes, we would like to have a 
prescription which relates the exact £ to the mean correla- 
tion function calculated from the linear theory. One might 
have naively imagined that £(a,x) should be related to 
(,l(cl,x). But one can convince oneself that the relationship 
is likely to be nonlocal by the following analysis: 

Recall that, the conservation of pairs of particles, gives 



an exact equation satisfied by the correlation function ( Pee- 
bles 19$0|): 



(2) 



where v(t,x) denotes the mean relative velocity of pairs at 
separation x and epoch t. Using the mean correlation func- 
tion £ and a dimensionless pair velocity h(a,x) = — (v/ax), 
equation (0) can be written as 
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(3) 



This equation can be simplified by first introducing the vari- 
ables 



A = In a, 



X = Inx, 



D(X,A)=ln(l + 



(4) 



in t erms of which we have ( Nityananda and Padmanabhan 
1994|) | 



^- h (A,X)^ = 3h(A,X) 



(5) 



Introducing further a variable F = D + 3X, equation Q 
can be written in a remarkably simple form as 



dF dF 

M -HA,x)- = o 



(6) 



The characteristic curves to this equation - on which F is 
a constant - are determined by (dX/dA) = — h(X, A) which 
can be integrated if h is known. But note that the charecter- 
istics satisfy the condition 



F = 3X + D = In [a: (1 + £)] = constant 
or, equivalently, 
x 3 (l + = l 3 



(7) 



(8) 



where I is another length scale. When the evolution is linear 
at all the relevant scales, f < 1 and I ~ x. As clustering 
develops, £ increases and x becomes considerable smaller 
than I. It is clear that the behaviour of clustering at some 
scale x is determined by the original linear power spectrum 
at the scale I through the "flow of information" along the 
charesteristics. This suggests that we should actually try to 



express the true correlation function £(a, x) in terms of the 
linear correlation function £l(o,1) evaluated at a different 
point. 

Let us see how we can do this starting from the quasi- 
linear regime. Consider a region surrounding a density peak 
in the linear stage, around which we expect the clustering 
to take place. It is well known that density profile around 
this peak can be described by 



p(x) « pbg[l + $(x)] 



(9) 



Hence the initial mean density contrast scales with the initial 
shell radius I as Si(l) oc in the initial epoch, when linear 
theory is valid. This shell will expand to a maximum radius 
of Xmax oc I / Si oc I / (,l(1) ■ In scale-invariant, radial collapse, 
models each shell may be approximated as contributing with 
a effective radius which is propotional to x max - Taking the 
final effective radius x as proportional to x ma x, the final 
mean correlation function will be 

M I 3 

i QL {x) oc P oc - oc (ta/ fr (t))3 « &(0 8 (10) 

That is, the final correlation function £ql at a; is the cube of 
initial correlation function at I where I 3 oc x 3 ^ oc x 3 £ql(x). 
This is in the form demanded by equation (||) if £ S> 1. Note 
that we did not assume that the initial power spectrum is a 
power law to get this result. 

In case the initial power spectrum is a power law, with 
ft oc x~(™ +3 \ then we find that 

| Qi oc x -3<"+3)/("+4) (ii) 

[If the correlation function in linear theory has the powerlaw 
form £l oc x~ a then the process desribed above changes the 
index from a to 3a/(l + a). We shall comment more about 
this aspect at the end of the paper.]. For the power law case, 
the same result can be obtained by more explicit means. For 
example, in power law models the energy of spherical shell 
will scale with its radius as some power which we write as 
E oc x 2 ~ b . Since M oc xf, it follows that the maximum 
radius reached by the shell scales as x max oc (M/E) oc x\ +h . 
Taking the effective radius as x = x e ff oc x] +b , the final 
density scales as 



M x 3 



-3b -36/(1+6) 
OC X; (XX ' v ' 



(12) 



In this quasilinear regime, c; will scale like the density and 
we get |q£ oc x~ 3b ^ 1+b K The index b can be related to n by 
assuming the the evolution starts at a moment when linear 
theory is valid. The gravitational potential energy [or the 
kinetic energy] scales as £ oc x~ ( - n+1 ^ in the linear theory. 
This may be seen as follows: The power spectrum for velocity 
field, P v (fc) in the linear regime is related to that of density 
by P v oc P(k)/k 2 oc k n ~ 2 . Hence the contribution to v 2 
in each logarithmic scale in k-space is k 3 P v /2Tv 2 oc fc n+1 oc 
x~( n+1 . Similarly, the gravitational potential energy due to 
fluctuations is 



o V 



(13) 



So the total energy in the initial configuration scales as 
2+' n+1 ' allowing us to determine b = n + 3. This shows 
that the correlation function in the quasilinear regime is the 
one given by equation ([H]) . 
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£nl{o,,x) oc [£ L (a,l)] 3/2 



(19) 
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Figure 1. Plot of f(a, x) against £l(<i, Z) for CDM model. The 
slopes in the three different regimes are indicated. The data points 
are for three different redshifts [0.1,0.5 and 1.0] and are based on 
the simulations described in Padmanabhan et al. (1995). 



The case with power law initial spectrum has no in- 
trisic scale, HQ— 1. It follows that the evolution has to 
be self similar and £ can only depend on q = xa~ 2 ^ n+3 \ 
This allows us to determine the a dependence of £ql by 
substituting q for x in equation (|ll|). We find 



(14) 



(15) 



i 0i (o, a! )oca 6/( ' ,+4) a: - s( " +8)/( " +4) 

Direct algebra shows that 

£,Q L (a,x) oc [^ L (a,l)f 

reconfirming the local dependence in a and nonlocal depen- 
dence in spatial coordinate. This result has no trace of orig- 
inal assumptions [spherical evolution, scale-invariant spec- 
trum ....] left in it and hence once would strongly suspect 
that it will have far general validity. 

Let us now proceed to the third and nonlinear regime. 
If we ignore the effect of mergers, then it seems reason- 
able that virialised systems should maintain their densities 
and sizes in proper coordinates, i.e. the clustering should 
be "stable". This would require the correlation function to 
have the form £,nl{cl,x) — a 3 F(ax). [The factor a 3 arising 
from the decrease in background density]. From our previ- 
ous analysis we expect this to be a function of (,l(o,, I) where 
I 3 » x 3 £ NL {a,x) 



Let us write this relation as 



(,Nh{a,x) = a F(ax) = U[£z,(a,l)] 



(16) 



where U[z] is an unknown function of its argument which 
needs to be determined. Since linear correlation function 
evolves as a 2 we know that we can write £l(o,,1) = a 2 Q[l 3 ] 
where Q is some known function of its argument. [We are 
using I rather than I in defining this function just for future 
convenience of notation]. In our case I 3 = x 3 ^nl(ci, x) = 
(ax) 3 F(ax) = r 3 F(r) where we have changed variables from 
(a, x) to (a, r) with r — ax. Equation [li] now reads 

a 3 F(r) = Ufoia, I)] = U[a 2 Q[l 3 ]] = U[a 2 Q[r 3 F(r)]] (17) 

Consider this relation as a function of a at constant r. 
Clearly we need to satisfy U[c\a 2 ] — c^a 3 where ci and C2 
are constants. Hence we must have 

U[z]ocz 3/2 . (18) 
Thus in the extreme nonlinear end we should have 



[Another way deriving this result is to note that if £ = 
a 3 F(ax), then h = 1. Integrating equation (^) with ap- 
propriate boundary condition leads to equation (^).] Once 
again we did not need to invoke the assumption that the 
spectrum is a power law. If it is a power law, then we get, 

3(n + 3) (2Q) 



£ N L(a,x) oc a { ""x 7 = 

(n + 5) 

This result is based on the assumption of "stable cluster- 



ing" and was originally derived by Peebles ( Peebles 1965). 
It can be directly verified that the right hand side of this 
equation can be expressed in terms of q alone, as we would 
have expected. 

Putting all our results together, we find that the non- 
linear mean correlation function can be expressed in terms 
of the linear mean correlation function by the relation: 

(for < 1, £ < 1) (21) 
(for 1< l L < 5.85, 1< f < 200) 
(for 5.85 < |i, 200 < f) 



fi(o,0 

£(a,x)= h(a,l) 3 
14.14fi,(a,Z) 



3/2 



The numerical coefficients have been determined by conti- 
nuity arguments. We have assumed the linear result to be 
valid upto £ = 1 and the virialisation to occur at £ ~ 200 
which is result arising from the spherical model. The ex- 
act values of the numerical coefficients can be obtained only 
from simulations. 

The true test of such a model, of course, is N-body sim- 
ulations and remarkably enough, simulations are very well 
represented by relations of the above form. Figure 1 shows 
the results of a CDM simulation ba sed on the investigations 
carried out in Padman abhan et.al (1995). This data ca n be 
fitted by the relations ( Bagla and Padmanabhan 1993| ) : 

£l(ci,1) (forfi < 1, f < 1) 

£(a, x) = h{a,lf (for 1< < 5, 1< f < 125) 

11.2£ L (a,l) 3/2 (for5<fi, 125 <0 (22) 

[The fact that numerical simulations show a correlation be- 
tween £(a, x) an d gi, (a, I) was originally pointed out by 
Hamilton et al. (1991) who, however, tried to give a mul- 
tiparameter fit to the data. This fit has somewhat obscured 
the simple physical interpretation of the result though has 
the virtue of being very accurate for numerical work.] 

A comparison of equations (|il]) and ( |22[ ) shows that 
the physical processes which operate at different scales are 
well represented by our model. In other words, the processes 
descibed in the quasilinear and nonlinear regimes for an in- 
dividual lump still models the average behaviour of the uni- 
verse in a statistical sense. It must be emphasised that the 
key point is the "flow of information" from I to x which is 
an exact result. Only when the results of the specific model 
are recast in terms of suitably chosen variables, we get a 
relation which is of general validity. It would have been, for 
example, incorrect to use spherical model to obtain relation 
between linear and nonlinear densities at the same location 
or to model the function h. With hindsight, it is clear why 
such attempts have not succeeded in the past. 

It may be noted that to obtain the result in the non- 
linear regime, we needed to invoke the assumption of sta- 
ble clustering which has not been deduced from any fun- 
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damental considerations. In case mergers of structures are 
impo rtant, one would consider this assumption to be sus- 



pect ( Padmanabhan, Cen, Ostriker and Summers 1995). We 
can, however, generalise the above argument in the following 
manner: If the virialised systems have reached stationarity 
in the statistical sense, the function h - which is the ratio 
between two velocities - should reach some constant value. 
In that case, one can integrate equation and obatin the 
result £nl = a ih F(a h x). A similar argument will now show 
that 



£jvz,(a, x) oc [£l(o, /)] 



3h/2 



(23) 



in the general case. For the power law spectra, one would 
get 



f (a, x) tx a' 3 1 ' ,h x 



;7 = 



3h(n + 3) 



(24) 



h{n + 3) 

Simulations are not accurate enough to fix the value of h; 
in particular, the asynptotic value of h could depend on n 
within the accuracy of the simulations. It may be possible 
to determine this dependence by modelling mergers in some 
simplified form. 

We conclude with two interesting speculations regard- 
ing the nonlinear stage. If h — 1 asymptotically, the correla- 
tion function in the extreme nonlinear end depends on the 
linear index n. One may feel that physics at highly nonlin- 
ear end should be independent of the linear spectral index 
n. This will be the case if the asymptotic value of h satisfies 
the scaling 



3c 
n + 3 



(25) 



in the nonlinear end with some constant c. Only high resso- 
lution numerical simulations can test this conjecture that 
h(n + 3) ^constant. 

Also note that the radial, scale invariant infall described 
in the quasilinear regime has the effect of changing the linear 
correlation function £l = a;~' n+3 - ) 
correlation function £ql = x 
what will be the effect of iterating this process N-times. It is 
easy to see that the index after N iterations can be expressed 
in the form: 



= x 
-36/(1+6) 



b to the quasilinear 
It is amusing to ask 



A N b 



IN 



An = 3 ; Bn = 



(26) 



1 + Bjvb' ' 2 

The fixed point, of course, is = 2 which is the only non- 
trivial fixed point for such an evolution [with the other, triv- 
ial, fixed point being zero]. If one could model the evolution 
as repeated application of this process, one would expect 
a continuum of scaling relations with the evolution being 
driven to a singular isothermal sphere. The quasilinear evo- 
lution doesnot change the x~ 2 profile, a r esult which was 
noted earlier in Bagla and Padmanabhan fll995j ). It is not 
clear whether the clustering can indeed be modelled using 
equation (p6"|). 

The relations obtained in this paper will, of course, have 
certain limitations on their validity. To begin with, we do 
expect a weak n-dependence in these relations due to av- 
eraging over peaks of different heights. This has been dis- 
cussed using a simple analy tic m odel, as well as numer ically , 
in T. Padmanabhan et.al, ( 1995 ) [Also see Mo et.al., (199J) 
for a similar discussion] . Secondly, the asymptotic behaviour 
will be sensitive to the value of Q. When f2 < 1, structures 



"freeze out" during the late stages of evolution and "stable 
clustering" is likely to be a reasonable assumption. Finally, 
models like HDM evolve in a manner very different from 
hierarchical models. In the former, small scale power is gen- 
erated by the breaking of long wavelength modes and the 
evolution is best modelled by instability of shell-like regions 
in the universe. Work is now in progress to generalise the 
ideas of the present paper for other models. 
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